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This  final  report  summarizes  our  major  accomplishments  on  the  research  under  AFOSR  Grant  FA9550- 
10-1-  0045  during  the  project  period  between  April  1,  2010  -  March  31,  2012. 

1  Background 

This  project  continues  our  previous  AFOSR  project  (Grant  No.  FA9550-08-1-0122)  to  extend  and  verify  our 
high  order  space-time  cell  vertex  scheme  (DG-CVS)  toward  solving  the  compressible  Navier-Stokes  equation. 

The  DG-CVS  method  integrates  the  best  features  of  the  space-time  Conservation  Element/Solution 
Element  (CE/SE)  method  [  ]  and  the  discontinuous  Galerkin  (DG)  method  [  ].  The  core  idea  is  to  construct  a 
staggered  space-time  mesh  through  alternate  cell-centered  CEs  and  vertex-centered  CEs  (cf.  Fig.  1  (right)) 
within  each  time  step.  Inside  each  SE  (cf.  Fig.  1  (left)),  the  solution  is  approximated  using  high-order 
space-time  DG  basis  polynomials.  The  space-time  flux  conservation  is  enforced  inside  each  CE  using  the 
DG  discretization.  The  solution  is  updated  successively  at  the  cell  level  and  at  the  vertex  level  within  each 
physical  time  step.  For  this  reason  and  the  method’s  DG  ingredient,  the  method  was  named  as  the  space-time 
discontinuous  Galerkin  cell- vertex  scheme  (DG-CVS). 

DG-CVS  equally  works  on  higher  dimensions  on  arbitrary  grids.  Figure  2  shows  the  conservation  elements 
and  solution  elements  on  quadrilateral  meshes  and  triangular  meshes.  Obviously,  the  definitions  of  CEs  and 
SEs  on  higher  dimensions  are  analogous  to  that  for  1-D  meshes  (cf.  Fig.  1).  Figure  3  demonstrates  the 
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Figure  1:  Solution  elements  (SEs)  and  conservation  elements  (CEs)  in  the  x  —  t  domain.  Left:  solution 
elements  and  right:  conservation  elements. 
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Figure  2:  Conservation  elements  (CEs)  and  solution  elements  (SEs)  in  the  x  —  y  —  t  domain.  First  row:  CEs 
for  rectangular  meshes,  second  row:  CEs  for  triangular  meshes,  third  row:  SEs  for  rectangular  meshes,  and 
fourth  row:  SEs  for  triangular  meshes. 
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Figure  3:  Dual  meshes  for  the  solution  updating  at  the  cell  level  (in  red)  and  the  vertex  level  (in  black). 
Left:  rectangular  mesh  and  right:  triangular  mesh. 


resulting  dual  mesh  at  the  cell  level  and  the  vertex  level  for  both  rectangular  meshes  and  triangular  meshes, 
respectively. 

The  main  features  of  DG-CVS  are  summarized  as  follows: 

•  based  on  space-time  formulation. 

•  arbitrary  high- order  accuracy  in  both  space  and  time.  Space  and  time  are  handled  in  a  unified  way 
based  on  space-time  flux  conservation  and  high-order  space-time  discontinuous  basis  functions.  This  is 
in  contrast  to  semi-discrete  methods  where  the  temporal  order  of  accuracy  is  limited  by  the  Backward 
Difference  Formula  (BDF)  or  the  multi-step  Runge-Kutta  method. 

•  Riemann- solver  free.  DG-CVS  does  not  need  a  (approximate)  Riemann  solver  to  provide  numerical 
fluxes  as  needed  in  finite  volume  or  traditional  DG  methods.  The  Riemann- solver- free  feature  offers 
two- fold  advantages.  First,  this  Riemann- solver- free  approach  eliminates  some  pathological  behaviors 
(e.g.,  carbuncle  phenomenon,  expansion  shocks,  etc.)  associated  with  some  Riemann  solvers.  Second, 
it  is  suitable  for  any  hyperbolic  PDE  systems  whose  eigenstructures  are  not  explicitly  known. 

•  reconstruction  free.  DG-CVS  solves  for  the  solution  and  its  all  spatial  and  temporal  derivatives  simul¬ 
taneously  at  each  space-time  node,  thus  eliminating  the  need  of  reconstruction. 

•  suitable  for  arbitrary  spatial  meshes.  The  CE  and  SE  definitions  in  DG-CVS  are  independent  of  the 
underlying  spatial  mesh  and  the  same  definitions  can  be  easily  extended  from  1-D  to  higher-dimensions 
(cf.  Fig.  2)  without  any  ambiguity.  Though  not  shown,  the  same  CE  and  SE  definitions  also  apply  to 
meshes  with  hanging  nodes. 

•  highly  compact  regardless  of  order  of  accuracy.  Only  information  at  the  immediate  neighboring  nodes 
will  be  needed  to  update  the  solution  at  the  new  time  level.  Compactness  eases  the  parallelization  of 
the  flow  solver. 

In  the  previous  project,  the  DG-CVS  method  has  been  developed  for  solving  hyperbolic  conservation  laws 
including  the  scalar  advection  equation  and  the  compressible  Euler  equations. 

2  Summary  of  Major  Accomplishments 

During  the  last  two  years,  we  have  accomplished  the  following: 

•  Verified  that  the  DG-CVS  method  solves  the  time- dependent  diffusion  equations  as  well  as  the  advection 
equations  without  any  special  care  on  diffusion  terms.  Thanks  to  the  staggered  space-time  conserva¬ 
tion  elements,  both  the  advective  and  diffusive  flux  are  continuous  and  unique  across  the  spatial  cell 
interface.  Therefore,  no  extra  reconstruction  or  recovery  or  ad  hoc  penalty  and  coupling  terms  are 
needed  to  avoid  the  “variational  crime”  and  ensure  the  consistency  of  the  variational  form  for  diffusive 
terms.  For  this  reason,  DG-CVS  is  conceptually  simper  than  other  existing  DG  methods  for  diffusion 
equations.  This  paves  the  way  to  solve  compressible  Navier-Stokes  equations. 
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•  Conducted  the  convergence  rate  for  diffusion  equations.  The  convergence  rate  is  L 2  optimal  when  the 
degree  of  basis  polynomials  is  odd  and  sub-optimal  when  the  degree  of  basis  polynomials  is  even.  In 
other  words,  DG-CVS  is  (p  +  l)th  order  accurate  for  the  solution  and  pth  order  accurate  for  the 
solution  gradients  for  odd  p.  When  p  is  even,  the  convergence  rate  of  DG-CVS  is  still  optimal  for 
advection  equations  but  sub-optimal  for  diffusion  equations.  In  practice,  one  can  choose  the  odd  p  for 
best  performance  when  solving  advection-diffusion  equations. 

•  Extended  to  more  practical  problems  involving  curved  boundaries.  The  curved  boundary  is  approxi¬ 
mated  by  high-order  geometric  polynomials  which  is  consistent  with  the  high-order  solution  polyno¬ 
mials.  This  is  indispensable  in  high-order  methods  to  avoid  spurious  entropy  generation  on  the  curved 
boundary. 

•  Extended  to  solve  the  shallow  water  equations.  The  DG-CVS  method  provides  an  alternative  Riemann- 
solver  free  approach  for  shallow  water  equations.  Being  able  to  solve  both  Euler  equations  and  shallow 
water  equations  confirms  that  DG-CVS  is  truly  a  working  Riemann- solver- free  approach  which  can  be 
used  to  solve  arbitrary  conservation  laws  where  eigenstructures  are  difficult  to  obtain  (e.g.,  magneto¬ 
hydrodynamics  equations)  in  the  future. 

•  Extended  to  solve  the  level  set  equation.  The  DG-CVS  method  is  able  to  resolve  the  interface  evolution 
governed  by  the  level  set  equation. 

•  Extended  to  solve  moving-mesh  problems.  Thanks  to  the  space-time  nature  of  the  DG-CVS  method, 
the  geometric  conservation  law  involving  moving  meshes  is  automatically  satisfied.  Therefore,  the 
DG-CVS  method  has  great  potentials  in  moving  boundary  problems  (e.g.,  fluid-structure  interaction). 

3  Description  of  the  Space-Time  Discontinuous  Galerkin  Cell- Vertex 
Scheme 

The  high-order  space-time  DG-CVS  method  is  developed  to  solve  the  equations  of  the  following  type: 

<9u  +  +  =  n  m 

dt  dx  dy  dx  dy  [  ’ 

where  u  is  the  conservative  state  vector,  f  and  g  are  advective  flux  vectors,  and  fv  and  gv  are  the  diffusive 

flux  vectors. 


3.1  Space-Time  Discontinuous  Galerkin  Formulation 

Following  the  idea  of  the  discontinuous  Galerkin  (DG)  method,  an  approximate  solution  is  sought  within 
each  space-time  solution  element  (SE),  denoted  as  K.  When  restricted  to  the  SE,  Uh  belongs  to  the  finite 
dimensional  space  U(K)  such  that 

Assume  all  the  components  of  the  conservative  variables  inside  each  SE  are  approximated  by  polynomials 
of  the  same  degree,  i.e 

N 

uh(x,t)-'52sj(l>j(x,t)  (2) 

3= 1 

where  {0j(x,  t)}^=1  are  some  type  of  space-time  polynomial  basis  functions  defined  within  the  solution 
element,  {s j}f=i  are  the  unknowns  to  be  determined  and  N  is  the  number  of  basis  functions  depending  on 
the  degree  of  the  polynomial  function. 

The  first  spatial  derivative  of  the  solution  can  be  expressed  in  terms  of  the  basis  functions  as  follows. 


duhhc,t)  d(f>j  duhht,t) 

irT¥Sj’  and^T“ 


N 


E 


(3) 


Currently,  the  polynomials  based  on  the  Taylor  expansions  are  used  as  the  basis  functions.  Note  that, 
for  the  Taylor  polynomials,  the  derivatives  of  the  polynomial  are  a  subset  of  the  original  Taylor  polynomials. 
This  allows  efficient  implementation  when  evaluating  integrals  involving  products  of  the  Taylor  polynomials. 
For  example,  the  quadratic  2-D  basis  functions  are  tabulated  in  Table  1. 
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Table  1:  Derivatives  of  quadratic  2-D  basis  functions. 
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Following  the  Galerkin  orthogonality  principle,  multiply  (1)  with  each  of  the  basis  functions  fa  (i  = 
1,  2,  •  •  •  ,  N)  and  integrate  over  a  space-time  CE  to  obtain  the  weak  form 


dg 

dy 


dfv 

dx 


dD  =  0, 


iml,  2,  ...,7V 


(4) 


where  Qk  is  the  space-time  conservation  element  (CE)  associated  with  the  solution  element  K.  Note  that 
the  conservation  element  is  identical  to  the  solution  element  except  for  the  volumeless  vertical  spike  in  the 
solution  element.  The  space-time  flux  conservation  in  weak  form  as  in  (4)  is  for  each  individual  space- 
time  conservation  element.  Therefore,  the  current  method  can  be  considered  as  a  space-time  discontinuous 
Galerkin  method. 

Integrating  (4)  by  parts  results  in 


/ 


dt  dx  [  v,+  dy[s  Sv[ 


L 


dQ  =  I  fatt-nd T 


(5) 


where  H  =  (fh  —  f^,  —  gj,  uh)  is  the  space-time  flux  tensor  and  n  =  (nx,  ny ,  nt)  is  the  outward  unit 

normal  of  the  CE  boundary,  i.e.  T  =  of  the  space-time  conservation  element  (CE)  under  consideration. 

Note  that  the  partial  integration  is  also  performed  on  the  time-dependent  term,  which  is  a  salient  difference 
between  space-time  DG  methods  and  semi-discrete  DG  methods.  As  can  be  seen,  the  formulation  in  (5) 
contains  both  the  volume  integral  and  the  surface  integral. 


3.2  Cell- Vertex  Solution  Updating  Strategy 

The  DG-CVS  inherits  the  core  idea  of  the  CE/SE  method  using  a  staggered  space-time  mesh  to  enforce 
the  space-time  flux  conservation.  However,  the  construction  of  the  staggered  space-time  slabs  in  DG-CVS 
deviates  from  the  CE/SE  method.  In  DG-CVS,  unknowns  are  stored  at  both  vertices  and  cell  centroids  of 
the  spatial  mesh,  and  the  solutions  at  vertices  and  cell  centroids  are  updated  at  different  time  levels  within 
each  time  step.  At  the  beginning  of  each  physical  time  step,  the  solution  is  assumed  known  at  the  vertices  of 
the  mesh,  either  given  as  the  initial  condition  or  obtained  from  the  previous  time  step.  Inside  each  new  time 
step,  the  solution  is  updated  in  two  successive  steps.  The  first  step  updates  the  solution  at  cell  centroids  at 
the  half-time  level  (tn+1/2)  based  on  the  known  vertex  solutions  at  the  previous  time  level  (£n).  The  second 
step  updates  the  solution  at  vertices  at  the  new  time  level  (£n+1)  based  on  the  known  cell  solutions  at  the 
previous  half-time  level  (tn+1/2).  The  same  process  is  repeated  for  new  time  steps. 

The  solution  updating  at  the  cell  level  or  the  vertex  level  is  based  on  the  key  equation  (5).  First  divide 
the  conservation  element  into  the  following  portions: 

•  the  interior  volume  Qk  where  the  solution  is  associated  with  the  space-time  node  at  the  new  time  level. 

•  the  top  surface  T^0p  where  the  solution  is  also  associated  with  the  space-time  node  at  the  new  time 
level. 
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Figure  4:  Illustration  of  space-time  flux  conservation  on  a  1-D  cell- level  CE.  Left:  stationary  mesh  and  right: 
moving  mesh. 


•  the  side  surfaces  Tg-qe  where  the  solution  is  associated  with  the  space-time  node  at  the  previous  time 
level. 

•  the  bottom  surfaces  Tp^  where  the  solution  is  also  associated  with  the  space-time  node  at  the  previous 
time  level. 


Correspondingly,  Eq.  (5)  can  be  rearranged  to  yield 
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(6) 


where  the  left  hand  side  contains  the  unknowns  and  the  right  hand  side  contains  the  known  information. 

Since  the  top  and  the  bottom  surface  of  the  CE  are  horizontal,  which  leads  to  •  n  =  uh  on  the  top 
face  and  Hh  •  n  =  —uh  on  the  bottom  face,  Eq.  (6)  can  be  simplified  to 


Qk  [  dt  dx 

f  <f>iUh6T 
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To  further  illustrate  the  idea  of  enforcing  the  space-time  flux  conservation,  consider  the  1-D  case  shown 
in  Fig.  4.  Suppose  the  solution  at  the  spacetime  node  (m  +  \,n  +  |)  is  to  be  determined.  Here,  m  and  n 
represents  the  spatial  and  the  temporal  locations  of  the  space-time  node,  respectively.  The  boundaries  of 
the  CE  associated  with  the  spacetime  node  (m  +  \,n  +  \)  is  divided  into  five  sections  T i,  T2,  T 3,  T4  and 
T5,  as  shown  in  Figure  4.  Among  these  sections,  T 1  belongs  to  the  SE  associated  with  node  (m  +  \,n  +  |) 
whose  solutions  are  to  be  determined,  T2  and  T3  the  SE  associated  with  node  (m,n)  and  T4  and  T5  the  SE 
associated  with  node  (m  +  l,n).  The  interior  volume  of  the  conservation  element  is  also  associated  with 
node  (m  +  |,n+  \).  Exactly  the  same  idea  can  be  applied  to  the  multidimensional  cases. 

As  can  be  seen,  with  this  staggered  space-time  cell-vertex  solution  updating  strategy,  no  Riemann- solver- 
typed  flux  functions  are  needed  for  the  interface  flux.  There  is  no  “left  state”  and  “right  state”  when  evaluating 
inter-cell  fluxes  as  Riemann  solvers  do.  We  can  see  the  Riemann- solver  free  DG-CVS  method  perfectly 
captures  all  flow  features  (shock  waves,  contact  discontinuities,  etc.)  without  pathological  phenomena  such 
as  the  expansion  shock. 

DG-CVS  also  solves  moving  mesh  problems  (Fig.  4  right)  without  special  care.  The  Geometrical  Con¬ 
servation  Law  (GCL)  which  is  vital  on  moving  meshes  is  automatically  satisfied  by  the  space-time  DG-CVS. 
This  is  because  the  mesh  speed  is  automatically  accounted  for  in  the  outward  unit  normal  n  =  (nx,ny,  nt)  of 
the  space-time  faces  of  the  CE.  Note  that  nt  is  in  general  nonzero  for  moving  meshes.  Therefore,  DG-CVS 
is  also  suitable  for  problems  with  moving  boundaries  and/or  r- typed  mesh  adaptation. 
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3.3  Solving  the  Local  Nonlinear  Equation  System 

Eq.  (7)  is  a  nonlinear  equation  system  for  each  space-time  node.  The  Newton-Raphson  iterative  method  is 
employed  to  solve  the  system.  To  apply  the  Newton-Raphson  method,  Eq.  (7)  is  first  written  in  the  form 

G(s)  =  0  (8) 

by  moving  the  right  hand  side  of  (7)  to  the  left  hand  side.  According  to  the  Newton-Raphson  method,  the 
solution  can  be  updated  iteratively  by  means  of 

s(r+l)  =s(t)  _  j-iG(s<t))  (9) 

where  the  superscript  r  represents  the  iteration  step.  The  Jacobian  matrix  J  = 

To  update  s  using  Eq.  (9),  one  can  first  compute  the  solution  increment  Js  by  solving 

JJs  =  -G(s(r))  (10) 

which  is  a  linear  equation  system  and  can  be  solved  by  the  LU-decomposition  method.  The  solution  is  then 
updated  via  s(r+1)  =  +  5s. 

Inside  each  Newton  iteration,  a  linear  equation  system  must  be  solved.  To  improve  the  efficiency,  the 
Jacobian  can  be  fixed  during  the  iteration  process.  In  the  current  implementation,  the  Jacobian  is  computed 
only  once  using  the  initial  solution  before  the  iteration  process.  The  Jacobian  is  then  LU-decomposed 
and  remain  unchanged  during  the  iteration.  Inside  each  iteration,  only  the  right  hand  side  of  (10)  is  updated 
and  only  the  back  substitution  operations  are  performed.  In  practice,  this  fixed  Jacobian  Newton  method 
works  very  well  and  only  slightly  increases  the  number  of  iterations.  However,  the  overall  efficiency  is  greatly 
improved. 

3.4  Quadrature-free  Implementation 

As  seen  in  Eq.  (7),  the  DG-CVS  formulation  involves  both  surface  integrals  and  volume  integrals.  An 
appropriate  quadrature  rule  (e.g.,  Gaussian  quadrature  rule)  is  typically  used  to  numerically  integrate  the 
surface  and  volume  integrals.  However,  quadrature  rules  are  expensive.  A  sufficiently  large  number  of 
quadratures  points  must  be  used  to  ensure  the  accuracy.  When  the  employed  basis  polynomial  is  of  degree 
p,  the  quadrature  rule  must  be  exact  for  polynomials  of  degree  2p  + 1  for  surface  integrals  and  must  be  exact 
for  polynomials  of  degree  2 p  for  volume  integrals  [3].  Therefore,  the  number  of  quadrature  points  increases 
rapidly  with  the  degree  of  polynomials  in  higher  dimensions.  Quadrature  free  implementation  is  desirable 
to  improve  the  efficiency  of  high  order  methods. 

To  allow  the  quadrature-free  implementation,  the  spatial  flux  vectors  must  be  expanded  in  terms  of  the 
basis  polynomials  similar  to  those  used  to  expand  the  conservative  variables.  For  example,  inviscid  fluxes 
can  be  expressed  as 

N  N 

fh  =  E  si^i’  and  8h  =  E  si  ‘A?  •  (n) 

j= i  j= i 

where  N  is  the  number  of  basis  polynomials  of  one  degree  higher  than  those  to  expand  the  conservative 
variables  as  in  Eq.  (2).  The  requirement  of  using  basis  polynomials  of  degree  p  +  1  is  necessary  to  ensure 
the  accuracy  of  the  scheme  in  the  case  of  nonlinear  fluxes.  The  method  based  on  the  Cauchy- Kovalewski 
(CK)  procedure  [4,  5]  is  especially  suitable  for  our  purpose.  In  the  current  implementation  of  DG-CVS,  the 
conservative  variables  are  expanded  using  the  Taylor  polynomials.  With  the  Taylor  polynomials,  the  spatial 
and  temporal  derivatives  of  the  solution  are  explicitly  available.  Using  the  CK  procedure,  the  space-time 
derivatives  of  the  flux  vectors  can  be  obtained  using  the  space-time  derivatives  of  the  conservative  variables. 
Taylor  polynomials  of  degree  p  have  the  following  general  form 

4>j  =  cxlymtn  (12) 

where  c  is  a  constant,  x  =  (x  —  xq)/5x,  y  =  (y  —  yo)/5y ,  and  t  =  (t  —  to) /St.  Here  (xo,  yo ,  to)  is  the  location 
of  the  space-time  node  whose  solutions  are  to  be  solved.  5x,  5y  and  5t  are  the  characteristic  sizes  of  the  local 
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space-time  conservation  element.  The  scaling  using  5x,  Sy  and  St  helps  to  improve  the  condition  number  of 
the  local  system.  In  (12),  /,  m,  and  n  are  integer  exponents  ranging  from  0  to  p  and  0  <  /  +  ra  +  n  <  p. 

Note  that,  for  the  Taylor  polynomials,  the  derivatives  of  the  polynomial  are  a  subset  of  the  original 
Taylor  polynomials.  Therefore,  the  products  and  all  have  the  forms  similar  to 

(12).  This  fact  allows  efficient  implementation  when  evaluating  integrals  involving  various  products  of  the 
Taylor  polynomials. 

Such  quadrature- free  integration  has  been  implemented  in  our  DG-CVS  solver  to  avoid  the  volume  integral 
and  the  integral  over  the  top  surface  of  the  CE  on  the  left  hand  side  of  Eq.  (7).  The  main  idea  is  to  use  the 
divergence  theorem  [6]  to  reduce  the  volume  integral  to  surface  integrals  and  further  to  line  integrals.  And 
finally  the  analytical  formulae  from  [7]  is  used  to  evaluate  the  line  integrals.  The  efficiency  of  the  integration 
is  significantly  increased  with  such  quadrature-free  technique.  The  benefit  is  higher  for  higher  p. 


3.4.1  Cauchy- Kovalewski  Procedure 

To  implement  quadrature-free  integration,  one  must  express  the  nonlinear  flux  vectors  in  terms  of  expan¬ 
sions  of  basis  polynomials  used  by  the  state  vector.  The  method  based  on  the  Cauchy- Kovalewski  (CK) 
procedure[5,  4]  is  especially  suitable  for  our  purpose.  In  the  current  implementation  of  DG-CVS,  the  con¬ 
servative  variables  are  expanded  using  the  Taylor  polynomials.  With  the  Taylor  polynomials,  the  spatial 
and  temporal  derivatives  of  the  solution  are  explicitly  available.  Using  the  CK  procedure,  the  space-time 
derivatives  of  the  flux  vectors  can  be  obtained  using  the  space-time  derivatives  of  the  conservative  variables. 

To  accomplish  this  goal,  several  methods  can  be  used.  The  first  possible  method  is  based  on  the  L 2 
projection[  ].  For  example,  to  obtain  the  expansion  of  pu2 ,  one  may  express  pu2  =  (pu^pu^  which  leads  to 
p(pu2)  =  (pu)(pu).  The  I/2  projection  method  states  that 

r  N  r 

/  (sPjU  fij')  =  c, bi(pu)(pu)dfl ,  *  =  1,2, •••,7V  (13) 

j=i  JQk 


(  2  A 

K !  . 


.  However,  the  left  hand  side  matrix  is  not 


which  is  an  N  x  N  linear  system  and  can  be  solved  for 

a  stationary  mass  matrix  since  p  appears  in  the  left  hand  side  of  Eq.  (13).  During  the  iterative  process,  the 
matrix  must  be  recomputed  when  p  is  updated.  Therefore,  the  projection  method  is  not  quite  efficient. 

Another  method  is  to  reconstruct  the  flux  vectors  using  the  flux  values  at  a  nodal  set  [9]  which  is  large 
enough  to  support  the  reconstruction  of  polynomials  of  degree  p  +  1.  The  flux  values  at  each  is  computed 
using  the  conservative  variables  at  that  node.  This  method  is  especially  suitable  for  elements  where  the 
Lagrangian  basis  functions  are  clearly  defined.  With  the  Lagrangian  polynomials,  the  functions  can  be 
expanded  using  the  nodal  values. 

Yet  another  method  based  on  the  Cauchy-Kovalewski  (CK)  procedure[5,  ]  is  especially  suitable  for  our 
purpose.  In  DG-CVS,  the  conservative  variables  are  expanded  using  the  Taylor  polynomials.  With  the 
Taylor  polynomials,  the  spatial  and  temporal  derivatives  of  the  solution  are  explicitly  available.  Using  the 
CK  procedure,  the  space-time  derivatives  of  the  flux  vectors  can  be  obtained  using  the  space-time  derivatives 
of  the  conservative  variables. 

Following  Dyson  ['  ]  and  Dumbser  et  al.  [1],  the  Cauchy-Kovalewski  procedure  which  is  based  on  the 
multidimensional  extension  of  the  Leibniz  rule  can  be  formulated  as 


r(a,b,c)if  da+6+Vi(^y,*)/2Qr,:M) 

C  dxaybtc 


(14) 


EEE 

i= 0  j= 0  k= 0  L 


C  \  3(a-d  +  (&-i)  +  (c-fc)  fa  Qi+3+k 

k  J  dxa~ldyb-i dtc~k  dxldyWtk 


where 


etc.  represents  the  binomial  coefficients. 


Note  that  Eq.  (14)  is  used  to  evaluate  the  space-time  derivatives  of  the  product  of  any  two  space- 
time  functions  /i  and  /2.  For  example,  when  /i  =  pu  and  /2  =  u ,  the  space-time  derivatives  of  pu 2  can 
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be  obtained  using  Eq.  (14).  However,  the  space-time  derivatives  of  u  is  not  available  since  u  is  not  a 
conservative  variable.  Therefore,  using  Eq.  (14)  alone  is  not  sufficient  for  our  goal.  Fortunately,  Dumbser  et 
al.  [4]  present  the  following  modified  Leibniz  rule  to  obtain  the  space-time  derivatives  of  from  the  known 
space-time  derivatives  of  f\  and  /1/2: 


1 

h 


r{a,b,c)rt  t  t  N  _  da+b+cf2(x,y,t) 
L**  (Ji/2,/1,/2)  -  - Q^ybjc - 


' h(x,y,t)f2{x,y,t)  _  £(a,b,c) 
dxaybtc 


U1J2) 


(15) 


where  c[a’b'c^ (/i,  fa)  is  evaluated  using  Eq.  (14)  except  that  the  last  term  of  the  sum  (i.e.  i  =  a,  j  =  b  and 
k  =  c)  is  set  to  be  zero,  namely, 


a  b  c 


EEE 


a 

i 


A-'* 


(A  ,/2) 


ga +b+c/1(X;y^)f2(X;y^) 

dxaybtc 


b  \  (  c  \  9(a-d+(&-i)+(c-fc)y1  ai+i+fc/2 

j  )  \  k  J  dxa~ldyb~^ dtc~k  dxldyidtk 


with 


Qa+b+cfe 


dxaybt 


c 


=  0 


(16) 


Therefore,  to  obtain  the  space-time  derivatives  of  u ,  we  just  need  to  let  /i  =  p,  =  u  and  /1/2  =  pu 
and  utilize  Eqs.  (15)  and  (16). 

In  this  paper,  Eqs.  (14-16)  are  used  to  construct  the  flux  expansions.  Based  on  the  above  formulations, 
the  Cauchy- Kovalewski  procedure  can  be  summarized  for  the  flux  vectors  of  the  Euler  equations: 


•  Step  1:  use  Eqs.  (15)  and  (16)  to  obtain  the  space-time  derivatives  of  u  and  v  from  the  known 
space-time  derivatives  of  p,  pu  and  pv. 

•  Step  2:  use  Eq.  (14)  to  obtain  the  space-time  derivatives  of  pu 2,  puv:  and  pv 2  from  the  known 
derivatives  of  pu,  pv,  u  and  v. 

•  Step  3:  the  space-time  derivatives  of  the  pressure  P  can  be  obtained  readily  since 


P  =  (7  —  1)  [pE  —  0.5(/m2  +  pv2)] 


where  pE  is  a  conservative  variable. 

•  Step  4:  the  space-time  derivatives  of  pH  can  be  readily  obtained  via  pH  =  pE  +  P. 

•  Step  5:  use  Eq.  (14)  to  obtain  the  space-time  derivatives  of  pHu  and  pHv  from  the  known  derivatives 
of  pH,  u  and  v. 

After  all  these  five  steps,  all  quantities  in  the  flux  vectors  can  be  expanded  in  terms  of  the  Taylor  expansions. 


3.4.2  Converting  volume  integrals  to  surface  integrals 

The  integrand  in  the  volume  integral  (cf.  (7))  has  the  following  general  form 

q(x,y,t)  =  qi(x,y)AxlAymAtn  (17) 

which  results  from  the  products  between  Taylor  polynomials.  Here,  qi(x,y)  is  a  non-polynomial  spatial 
function.  For  instance,  in  some  applications,  the  source  term  cannot  be  explicitly  expressed  as  an  polynomial 
function.  So  q\{x,y)  is  kept  here  without  loss  of  generality.  In  (17),  l,  m,  and  n  are  integer  exponents.  For 

notational  simplicity,  denote  x  =  Ax  =  x  —  xq,  y  —  Ay  =  y  —  y$,  and  t  =  At  =  t  —  to-  Then  (17)  becomes 

q(x,y,t )  =  qi(x,y)xlymin  (18) 

Computing  volume  integrals  with  integrands  of  the  form  as  in  (18)  is  equivalent  to  computing  the  moments 
of  arbitrary  order  of  polyhedra.  In  the  current  case,  the  polyhedra  are  polygonal  cylinders.  The  method  of 
using  the  divergence  theorem  explained  in  [i  ]  is  adopted  to  reduce  the  volume  integral  to  surface  integrals. 
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To  evaluate  a  scalar  integral  inside  a  space-time  volume  fl,  i.e.  qdfl,  we  first  construct  a  vector  function 
Q (x,y,t)  =  Q i ii  +  Q 2^2  +  Q3I3  where  ii,  i 2  and  i3  denote  the  unit  vectors  along  the  x,  y ,  and  t-directions, 
respectively,  such  that 

q  =  V  •  Q  (19) 

With  the  aid  of  the  auxiliary  function  Q,  one  can  convert  the  volume  integral  into  a  surface  integral  via 
the  divergence  theorem  as  follows. 


/  q(x,y,t)dTl=  /  V  •  Qdfl  =  /  Q  •  ndT 

J  Q  J  Q  J  dQ 

Without  loss  of  generality,  we  can  define  the  following  for  Q 

Qi  =  0,  Q2  =  0,  and  Q3(x,  y,t)  =  J  q(x,  y,  t)dt  +  c(x,  y) 


(20) 


(21) 


to  satisfy  (19).  In  (21),  f  qdt  is  an  indefinite  integral  with  respect  to  t  and  c(x,y)  is  an  arbitrary  function 
independent  of  t. 

Substituting  (21)  into  (20)  to  obtain 

/  q(x,y,t)dQ  =  /  Q3ntdA  =  I  (  q(x,  y,  t)dt )  ntdA  +  /  c(x,y)ntdA  (22) 

Jn  Jon  Jan  \J  J  Jon 

where  nt  is  the  t-component  of  the  outward  unit  normal  of  the  surface. 

fdn  c(x,  y)nt dT  equals  zero  for  enclosed  surfaces  since  c(x,  y)  is  independent  of  t  and  results  in  zero  net 
contributions.  Therefore,  it  suffices  to  evaluate  fQ  qdfl  via 

[  q(x,y,t)  dQ=  [  (  [  q(x,y,t)dt)  ntdA  (23) 

Jn  Jan  \J  J 

Since  nt  =  1  and  nt  =  —  1  on  the  top  and  bottom  faces  of  the  CE,  respectively,  (23)  can  be  computed  as 

J  q(x,  y,  t)dQ  =  J  q(x,y,t)dtj  dA 

+  J  ( J  q(x,y,t)dtjntdA-  J  ( f  q(x,y,t)dtj  dA  (24) 

rside  v  /  rbott  v  ' 

Furthermore,  on  the  top  surface,  t  =  At  =  0,  therefore,  the  first  integral  on  the  right  hand  side  of  Eq. 
(24)  is  zero.  (24)  reduces  to 

J^q(x,y,t)dn  =  J  q(x,y,t)dtSjntdA-  q(x,y,t)dtj  dA  (25) 

side  bott 

For  simple  polynomials,  the  indefinite  integral  can  be  evaluated  analytically,  i.e. 

J  q(x,y,t)dt  =  qi(x,y)xlym  J  tndt  =  j-<?i(a;,  y)xlymtn+1  (26) 


Substituting  (26)  into  (25)  and  considering  the  fact  that  t  is  constant  on  the  horizontal  bottom  face  of 
the  CE  leads  to 

[  q(x,y,t)dfl  =  — J—  [  (q1(x,y)xlymtn+1)  ntdA 

Jn  n+Urside 

-  — j-r  tn+1  [  qi(x,y)xlymdA  (27) 

n  +  1  Jrbott 

As  can  be  seen,  the  volume  integral  q(x ,  y ,  t)dQ  has  been  reduced  to  a  surface  integral  (27)  which  can 

be  integrated  numerically  using  the  Gaussian  rule. 

The  integrals  in  (27)  do  not  need  to  be  computed  in  an  exclusive  subroutine.  Instead,  they  can  be 
computed  in  the  same  subroutine  where  the  surface  integrals  on  the  right  hand  side  of  Eq.  (7)  are  evaluated 
using  the  Gaussian  quadrature  rule.  Therefore,  the  extra  cost  of  computing  (27)  is  trivial. 
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4  Major  Progress  #1:  Solving  Advection-Diffusion  Equations 


4.1  Grid  Convergence  Study 

An  important  measurement  of  the  performance  of  any  high-order  method  is  its  convergence  order  of  accuracy. 
In  the  context  of  finite  element-based  methods,  when  basis  polynomials  of  degree  p  is  used,  and  if  the 
solution  u  and  the  solution  gradient  ux  are  approximated  using  Eqs.  (2)  and  (3),  respectively,  then  the 
optimal  convergence  orders  are  p  +  1  for  the  solution  u  and  p  for  the  solution  gradient.  Theoretical  analysis 
of  convergence  orders  of  high-order  methods  is  difficult,  it  not  impossible.  However,  the  convergence  order 
can  be  easily  determined  numerically  via  the  grid  convergence  study  with  appropriate  error  norms.  In  this 
sub-section,  the  grid  convergence  study  is  conducted  on  the  1-D  advection-diffusion  equation  and  the  2-D 
heat  equation.  The  convergence  behavior  of  DG-CVS  for  the  advection  equation  and  the  diffusion  equation 
will  be  compared. 

To  determine  the  numerical  convergence  order,  the  following  error  norms  are  defined: 


^oo  (e) 

h(e) 


max|uft(Xi)  -uexact(Xi) | 


N 


—  —  ^exact(x,))2 


i= 1 


i2(c)  = 


(28) 


where  nv  is  the  number  of  vertices  in  the  computational  domain,  uh(~Xi)  is  the  computed  numerical  solution 
at  ith  vertex  and  uexac^(jCi)  is  the  analytical  solution  at  the  ith  vertex,  fli  is  the  spatial  domain  associated 
with  the  ith  vertex  and  \Q\  is  the  size  of  the  entire  computational  domain.  Here,  and  l2  norms  are 
evaluated  at  the  discrete  location  of  vertices  and  L2  norm  is  obtained  by  integrating  the  continuous  solution 
within  the  spatial  domain  associated  with  each  vertex. 


4.1.1  Convergence  Orders  on  1-D  Meshes 

The  first  test  is  to  solve  the  following  1-D  linear  scalar  advection-diffusion  equation 


du  du  d2u 


dt  +  a  dx  dx 2 
u(x,  0)  =  Uo(x) 


0,  -1<x<1 

sin(7nr),  periodic  b.c. 


The  analytical  solution  is  given  as 

^exact  =  e  sin(7r(x  —  at)) 


(29) 


The  following  two  cases  are  tested: 

•  pure  advection  equation,  a  =  1.0,  v  =  0. 

•  heat  equation,  a  =  0,  v  m  1.0. 

For  each  case,  four  meshes  (number  of  cells  nc  =  10,  20,  40,  and  80)  are  used  for  varying  degrees  of  basis 
polynomials  (p  =  1,  2,  3,  and  4). 

For  the  advection  case,  the  time  step  is  chosen  as  St  =  crSx/a  where  Sx  is  the  cell  interval  and  the  Courant 
number  a  =  0.5,  0.3125,  0.25  and  0.25  for  pi,  p2,  p3  and  p4  cases,  respectively.  For  the  stability  limits  in 
the  case  of  advection  equations,  one  can  refer  to  our  earlier  paper[l  ].  All  cases  are  computed  up  to  t  =  1.0. 
The  numerical  convergence  orders  using  the  Z^,  l2  and  L2  norms  are  recorded  in  Tables  2,  3  and  4.  All  three 
tables  show  that  the  numerical  convergence  orders  are  p  +  1  for  u  and  p  for  ux  for  all  p’s.  The  convergence 
orders  are  optimal. 

For  the  diffusion  case,  the  time  step  is  chosen  as  St  =  0.1  Sx1  /v  for  all  pi  —  p4  cases.  All  cases  are 
computed  up  to  t  =  0.1.  The  numerical  convergence  orders  based  on  various  norms  are  recorded  in  Tables  5 
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Table  2:  1-D  advection  equation,  a  =  1.0,  t  =  1.0.  Numerical  convergence  order  determined  by  norm. 


nc  10 

nc  = : 

20 

nc  =  ■ 

40 

nc  —  ; 

80 

p 

variable 

loo  (D 

loo  (e) 

order 

loo  (e) 

order 

loo  (e) 

order 

comments 

1 

u 

2.11E-02 

4.88E-03 

2.113 

1.13E-03 

2.112 

2.70E-04 

2.064 

optimal 

Ux 

1.08E-01 

4.38E-02 

1.301 

2.09E-02 

1.071 

1.03E-02 

1.016 

optimal 

2 

u 

2.13E-04 

2.77E-05 

2.946 

3.51E-06 

2.979 

4.41E-07 

2.995 

optimal 

ux 

2.89E-02 

7.18E-03 

2.008 

1.79E-03 

2.001 

4.48E-04 

2.000 

optimal 

3 

u 

2.56E-05 

1.58E-06 

4.019 

9.71E-08 

4.023 

6.04E-09 

4.007 

optimal 

Ux 

4.04E-04 

6.03E-05 

2.742 

7.80E-06 

2.951 

9.83E-07 

2.989 

optimal 

4 

U 

4.18E-07 

1.28E-08 

5.031 

3.90E-10 

5.037 

1.25E-11 

4.965 

optimal 

Ux 

5.59E-05 

3.51E-06 

3.992 

2.21E-07 

3.991 

1.39E-08 

3.990 

optimal 

Table  3:  1-D  advection  equation,  a  =  1.0,  t  =  1.0.  Numerical  convergence  order  determined  by  I2 


nc  =  10 

nc  = : 

20 

nc  1=  - 

40 

nc  =  ■ 

80 

P 

variable 

/2(e) 

order 

h(e) 

order 

order 

comments 

1 

u 

1.48E-02 

3.36E-03 

2.134 

7.88E-04 

2.095 

1.90E-04 

2.055 

optimal 

Ux 

7.71E-02 

3.07E-02 

1.330 

1.46E-02 

1.071 

7.25E-03 

1.010 

optimal 

2 

u 

1.61E-04 

2.02E-05 

2.994 

2.52E-06 

3.006 

3.14E-07 

3.006 

optimal 

Ux 

2.13E-02 

5.20E-03 

2.037 

1.28E-03 

2.018 

3.19E-04 

2.009 

optimal 

3 

u 

1.76E-05 

1.09E-06 

4.013 

6.79E-08 

4.009 

4.24E-09 

3.999 

optimal 

Ux 

2.86E-04 

4.18E-05 

2.774 

5.45E-06 

2.940 

6.91E-07 

2.981 

optimal 

4 

u 

3.12E-07 

9.28E-09 

5.069 

2.79E-10 

5.057 

8.93E-12 

4.964 

optimal 

Ux 

4.13E-05 

2.54E-06 

4.022 

1.58E-07 

4.008 

9.89E-09 

3.999 

optimal 

Table  4:  1-D  advection  equation,  a  =  1.0,  t  =  1.0.  Numerical  convergence  order  determined  by  L2 


nc  —  10 

nc  = : 

20 

nc  =  - 

40 

nc  =  1 

80 

P 

variable 

L2(e) 

L2(e) 

order 

L2(e) 

order 

L2(e) 

order 

comments 

1 

u 

1.17E-02 

2.79E-03 

2.064 

6.87E-04 

2.020 

1.71E-04 

2.006 

optimal 

Ux 

4.04E-01 

2.03E-01 

0.993 

1.02E-01 

0.998 

5.09E-02 

0.999 

optimal 

2 

u 

5.73E-04 

7.21E-05 

2.990 

9.03E-06 

2.997 

1.13E-06 

2.999 

optimal 

ux 

3.66E-02 

9.21E-03 

1.990 

2.31E-03 

1.998 

5.77E-04 

1.999 

optimal 

3 

u 

2.55E-05 

1.64E-06 

3.961 

1.03E-07 

3.988 

6.47E-09 

3.997 

optimal 

ux 

2.31E-03 

2.94E-04 

2.975 

3.69E-05 

2.992 

4.62E-06 

2.998 

optimal 

4 

u 

8.02E-07 

2.50E-08 

5.005 

7.80E-10 

5.002 

2.46E-11 

4.986 

optimal 

Ux 

1.00E-04 

6.25E-06 

3.999 

3.90E-07 

4.002 

2.45E-08 

3.996 

optimal 

norm. 


norm. 
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Table  5:  1-D  heat  equation,  v  —  1.0,  t  =  0.1.  Numerical  convergence  order  determined  by  norm. 


nc  =  10 

nc  = : 

20 

nc  =  40 

nc  =  1 

80 

V 

variable 

Joo(e) 

^00  (e) 

order 

loo(e) 

order 

Joo(e) 

order 

comments 

1 

u 

9.72E-03 

2.62E-03 

1.894 

6.58E-04 

1.991 

1.65E-04 

1.998 

optimal 

ux 

3.04E-02 

7.82E-03 

1.959 

1.97E-03 

1.990 

4.93E-04 

1.997 

super-optimal 

2 

u 

1.32E-03 

3.70E-04 

1.838 

9.40E-05 

1.977 

2.36E-05 

1.994 

sub-optimal 

Ux 

3.17E-02 

8.28E-03 

1.936 

2.09E-03 

1.983 

5.25E-04 

1.996 

optimal 

3 

u 

7.29E-05 

5.10E-06 

3.836 

3.24E-07 

3.976 

2.03E-08 

3.994 

optimal 

Ux 

5.08E-05 

3.66E-06 

3.792 

2.37E-07 

3.948 

1.50E-08 

3.987 

super-optimal 

4 

u 

5.14E-06 

3.23E-07 

3.993 

1.99E-08 

4.020 

1.24E-09 

4.006 

sub-optimal 

Ux 

2.79E-05 

1.77E-06 

3.974 

1.11E-07 

3.992 

6.98E-09 

3.998 

optimal 

Table  6:  1-D  heat  equation,  v  —  1.0,  t  =  0.1.  Numerical  convergence  order  determined  by  I2  norm. 


nc  =  10 

nc  =  : 

20 

nc  =  40 

nc  =  i 

80 

V 

variable 

h(e) 

h(e) 

order 

12(e) 

order 

h(e) 

order 

comments 

1 

u 

6.89E-03 

1.81E-03 

1.933 

4.60E-04 

1.974 

1.16E-04 

1.989 

optimal 

Ux 

2.25E-02 

5.66E-03 

1.988 

1.41E-03 

2.006 

3.51E-04 

2.006 

super-optimal 

2 

u 

9.38E-04 

2.55E-04 

1.877 

6.57E-05 

1.959 

1.66E-05 

1.985 

sub-optimal 

Ux 

2.34E-02 

5.99E-03 

1.965 

1.50E-03 

2.000 

3.74E-04 

2.004 

optimal 

3 

u 

5.16E-05 

3.52E-06 

3.875 

2.26E-07 

3.959 

1.43E-08 

3.985 

optimal 

Ux 

3.75E-05 

2.65E-06 

3.822 

1.70E-07 

3.964 

1.07E-08 

3.995 

super-optimal 

4 

u 

3.65E-06 

2.23E-07 

4.032 

1.39E-08 

4.002 

8.71E-10 

3.997 

sub-optimal 

Ux 

2.06E-05 

1.28E-06 

4.003 

7.98E-08 

4.008 

4.97E-09 

4.006 

optimal 

and  6  and  7.  Comparing  with  the  optimal  convergence  rates  for  the  advection  equation,  one  can  see  that  the 
convergence  rates  for  the  diffusion  equation  appear  inconsistent.  For  clarity’s  sake,  the  following  observations 
from  Tables  5  and  6  and  7  are  summarized: 

•  When  p  is  odd, 

—  the  convergence  rates  based  on  the  /oo  and  I2  norms  are  optimal  for  u  and  super-optimal  for  ux. 

—  the  convergence  rates  based  on  the  L2  norm  are  optimal  for  both  u  and  ux. 

•  When  p  is  even,  the  convergence  rates  based  on  all  norms  are  sub-optimal  for  u  and  optimal  for  ux. 

Since  /-norms  do  not  produce  the  same  convergence  orders  for  the  heat  equation  as  those  for  the  advection 
equation  no  matter  whether  p  is  odd  or  even,  while  the  convergence  orders  based  on  the  Z/2-norm  are  optimal 
for  both  the  heat  equation  and  the  advection  equation  when  p  is  odd,  this  indicates  that  Z/2-norm  is  a  more 
appropriate  norm  for  determining  the  convergence  order.  We  will  focus  on  Z/2-norm  in  the  remaining  tests. 

The  inconsistent  convergence  behavior  between  odd  degree  and  even  degree  approximations  has  also  been 
reported  in  the  methods  of  many  other  researcher  [11,  12,  13,  14,  15].  It  is  also  interesting  to  note  that  the 
first  version  of  the  central  LDG  method  by  Liu  et  al.[l(  ]  is  sub-optimal  first  order  accurate  for  p  =  1  and 
optimal  ( p  +  l)th  order  accurate  for  p  >  1. 

4.1.2  Convergence  Orders  on  2-D  Structured  and  Unstructured  Meshes 

To  investigate  the  convergence  behavior  of  DG-CVS  for  2-D  diffusion  equations,  rectangular  and  unstructured 
triangular  meshes  with  various  resolutions  are  used.  Figure  5  shows  the  coarsest  rectangular  and  unstructured 
triangular  meshes  used  in  the  test.  The  coarsest  rectangular  mesh  is  composed  of  10  x  10  rectangular  cells 
and  is  designated  as  “qua- 10”.  The  coarsest  triangular  mesh  is  designated  as  “tri-10”  whose  edge  resolution 
is  comparable  to  that  of  mesh  qua- 10.  The  meshes  will  be  refined  isotropically  several  times  in  the  grid 
convergence  study,  resulting  in  a  series  of  meshes  designated  as  “qua-20”,  “qua-40”,  “qua-80”,  “tri-20”,  “tri-40” 
and  “tri-80”,  respectively. 
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Table  7:  1-D  heat  equation,  v  =  1.0,  t  =  0.1.  Numerical  convergence  order  determined  by  L2  norm. 


nc  =  10 

nc  = 

20 

nc  = 

40 

nc  = 

80 

V 

variable 

L2(e) 

L2(e) 

order 

L2(e) 

order 

L2(e) 

order 

comments 

1 

u 

4.90E-03 

1.24E-03 

1.982 

3.11E-04 

1.995 

7.78E-05 

1.999 

optimal 

ux 

1.49E-01 

7.50E-02 

0.994 

3.75E-02 

0.999 

1.88E-02 

1.000 

optimal 

2 

u 

1.22E-03 

2.78E-04 

2.130 

6.75E-05 

2.042 

1.68E-05 

2.011 

sub-optimal 

Ux 

1.54E-02 

3.93E-03 

1.965 

9.90E-04 

1.991 

2.48E-04 

1.998 

optimal 

3 

u 

3.70E-05 

2.45E-06 

3.918 

1.55E-07 

3.978 

9.76E-09 

3.994 

optimal 

ux 

6.50E-04 

8.09E-05 

3.006 

1.01E-05 

3.000 

1.26E-06 

3.000 

optimal 

4 

u 

3.92E-06 

2.30E-07 

4.094 

1.41E-08 

4.025 

8.77E-10 

4.007 

sub-optimal 

ux 

3.75E-05 

2.15E-06 

4.128 

1.31E-07 

4.034 

8.14E-09 

4.008 

optimal 

Figure  5:  Two-dimensional  meshes  used  in  the  tests.  Left:  rectangular  mesh.  Right:  unstructured  triangular 
mesh. 
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Table  8:  2-D  heat  equation  with  sinusoidal  initial  solution  on  rectangular  meshes,  v  =  1.0,  t  =  0.1.  Numerical 
convergence  order  determined  by  the  L 2  norm. 


V 

variable 

qua- 10 

qua-20 

qua-40 

qua- 80 

comments 

L2(e) 

L2(e) 

order 

L2(e) 

order 

L2(e) 

order 

1 

u 

ux  and  uy 

4.27e-03 

7.81e-02 

1.09e-03 

3.94e-02 

1.970 

0.987 

2.73e-04 

1.98e-02 

1.997 

0.993 

6.84e-05 

9.89e-03 

1.997 

1.001 

optimal 

optimal 

2 

u 

ux  and  uy 

8.37e-04 

1.24e-02 

2.02e-04 

3.21e-03 

2.051 

1.950 

5.00e-05 

8.11e-04 

2.014 

1.985 

1.25e-05 

2.03e-04 

2.000 

1.998 

sub-optimal 

optimal 

3 

u 

ux  and  uy 

5.43e-05 

1.28e-03 

3.73e-06 

1.58e-04 

3.864 

3.018 

2.39e-07 

1.97e-05 

3.964 

3.004 

1.51e-08 

2.47e-06 

3.984 

2.996 

optimal 

optimal 

4 

u 

ux  and  uy 

2.92e-06 

1.01e-04 

1.68e-07 

6.06e-06 

4.119 

4.059 

1.05e-08 

3.75e-07 

4.000 

4.014 

6.54e-10 

2.34e-08 

4.005 

4.002 

sub-optimal 

optimal 

Table  9:  2-D  heat  equation  with  sinusoidal  initial  solution  on  triangular  meshes,  v  ==  1.0,  t  ==  0.1.  Numerical 
convergence  order  determined  by  the  L2  norm. 


tri-10 

tri-20 

tri-40 

tri-80 

V 

variable 

L2(e) 

L2(e) 

order 

L2(e) 

order 

L2(e) 

order 

comments 

1 

u 

4.82e-03 

1.25e-03 

1.947 

3.22e-04 

1.957 

8.12e-05 

1.988 

optimal 

ux 

7.26e-02 

3.64e-02 

0.996 

1.82e-02 

1.000 

9.13e-03 

0.995 

optimal 

Uy 

7.25e-02 

3.64e-02 

0.994 

1.82e-02 

1.000 

9.13e-03 

0.995 

optimal 

2 

u 

1.37e-03 

3.46e-04 

1.985 

8.65e-05 

2.000 

2.17e-05 

1.995 

sub-optimal 

ux 

9.07e-03 

2.30e-03 

1.979 

5.80e-04 

1.988 

1.46e-04 

1.990 

optimal 

Uy 

9.03e-03 

2.29e-03 

1.979 

5.78e-04 

1.986 

1.45e-04 

1.995 

optimal 

3 

u 

3.14e-05 

2.02e-06 

3.958 

1.28e-07 

3.980 

8.07e-09 

3.987 

optimal 

ux 

7.30e-04 

9.23e-05 

2.983 

1.16e-05 

2.992 

1.46e-06 

2.990 

optimal 

Uy 

7.30e-04 

9.23e-05 

2.983 

1.16e-05 

2.992 

1.46e-06 

2.990 

optimal 

4 

U 

2.15e-06 

1.33e-07 

4.015 

8.35e-09 

3.994 

5.27e-10 

3.986 

sub-optimal 

^X 

4.62e-05 

2.96e-06 

3.964 

1.88e-07 

3.977 

1.19e-08 

3.982 

optimal 

Uy 

4.62e-05 

2.96e-06 

3.964 

1.88e-07 

3.977 

1.19e-08 

3.982 

optimal 

The  following  2-D  heat  equation  with  sinusoidal  initial  solution  is  solved  using  DG-CVS. 


du  ( d2u  d2u\ 

dt  \  dx 2  dy 2  ) 

u{x:y,  0) 


0,  —  1  <  x,  y  <  1 

sin(7r(x  +  2/)),  periodic  b.c. 


The  analytical  solution  is  given  as 

uexact  =  sin(7r(a;  +  y)) 


(30) 


The  time  step  is  chosen  to  be  St  =  ah2  jv  where  h  is  the  local  characteristic  size  of  the  element  and 
a  =  0.2  for  pi  cases  and  a  =  0.1  for  all  other  cases.  Table  8  and  Table  9  show  the  convergence  orders  on 
rectangular  meshes  and  triangular  meshes,  respectively.  Both  tables  show  the  same  convergence  rates  which 
indicates  the  convergence  orders  are  independent  on  the  type  of  the  spatial  mesh.  In  addition,  the  same 
convergence  rates  as  those  in  Table  7  for  the  1-D  heat  equation  are  observed,  that  is,  the  convergence  orders 
are  optimal  for  both  u  and  its  gradient  when  p  is  odd,  and  sub-optimal  for  u  and  optimal  for  ids  gradients 
when  p  is  even. 


4.2  More  Tests  on  2D  Scalar  Advection-Diffusion  Equation 

In  this  sub-section,  we  provide  two  more  test  cases.  Both  cases  involve  time  dependent  boundary  conditions. 
The  first  case  is  the  2-D  heat  equation  with  the  delta  initial  solution,  and  the  second  case  is  the  2-D  nonlinear 
viscous  Burgers  equation. 
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Figure  6:  Solution  of  2-D  heat  equation  with  the  delta  initial  solution.  Top  row:  solutions  on  the  rectangular 
mesh  qua-50  and  bottom  row:  solutions  on  the  triangular  mesh  tri-50.  Left  column:  solution  at  t  =  0.01, 
middle  column:  solution  at  t  =  0.05  and  right  column:  solution  at  t  =  0.2. 


From  the  previous  tests,  one  can  conclude  that  the  convergence  orders  based  on  the  L 2  norm  is  optimal 
for  both  the  solution  and  its  gradients  when  p  is  odd.  Besides,  this  optimality  holds  for  both  the  advection 
and  diffusion  equations.  Since  in  practice,  the  governing  equations  often  involve  both  advection  and  diffusion 
simultaneously,  it  justifies  to  employ  basis  polynomials  of  odd  degrees  for  best  and  consistent  convergence 
rates.  Therefore,  in  the  remaining  test  cases,  only  results  of  odd  p  will  be  presented. 


4.2.1  2-D  Heat  Equation  with  the  Delta  Initial  Solution 

This  case  is  to  solve  the  following  2-D  heat  equation  using  DG-CVS. 


du  f  d2u  d2u\ 
dt  dy2  ) 

u{x,y,  0) 


0,  —  1  <  x,  y  <  1 

Uq 


(31) 


where  uo  is  the  delta  function  at  the  origin  of  the  domain  (0,0).  The  solution  to  Eq.  (31)  with  such  an 
initial  condition  is  called  the  fundamental  solution  of  the  heat  equation[17].  The  analytical  solution  is  given 
as 


uexact  -4^ 


i±vi 


The  boundary  conditions  on  the  four  boundaries  are  time  varying  depending  on  the  above  analytical  solution 
formula. 

The  time  step  is  chosen  to  be  St  =  crh2  where  h  is  the  local  characteristic  size  of  the  element  and  a  =  0.2 
for  the  pi  case  and  a  =  0.1  for  the  p3  case.  Figure  6  shows  the  carpet  view  of  the  solution  at  three  instants, 
t  =  0.01,  t  =  0.05,  and  t  =  0.2,  on  both  the  rectangular  and  triangular  meshes,  respectively. 


16 


Table  10:  2-D  heat  equation  with  delta  initial  solution  on  rectangular  meshes,  t  =  0.05.  Numerical  conver¬ 
gence  order  determined  by  L 2  norm. 


p 

variable 

qua- 10 

qua-20 

qua-40 

qua-80 

comments 

L2(e) 

L2(e) 

order 

L2(e) 

order 

L2(e) 

order 

1 

u 

ux  and  uy 

1.20e-02 

2.53e-01 

3.12e-03 

1.28e-01 

1.943 

0.983 

7.88e-04 

6.43e-02 

1.985 

0.993 

1.98e-04 

3.22e-02 

1.993 

0.998 

optimal 

optimal 

3 

u 

ux  and  uy 

4.53e-04 

6.82e-03 

3.33e-05 

7.78e-04 

3.766 

3.132 

2.19e-06 

9.38e-05 

3.927 

3.052 

1.38e-07 

1.15e-05 

3.988 

3.028 

optimal 

optimal 

Table  11:  2-D  heat  equation  with  delta  initial  solution  on  triangular  meshes,  t  =  0.05.  Numerical  convergence 
order  determined  by  L2  norm. 


tri-10 

tri-20 

tri-40 

tri-80 

P 

variable 

L2(e) 

L2(e) 

order 

L2(e) 

order 

L2(e) 

order 

comments 

1 

u 

1.15e-02 

3.02e-03 

1.929 

7.71e-04 

1.970 

1.95e-04 

1.983 

optimal 

ux 

2.23e-01 

1.12e-01 

0.994 

5.64e-02 

0.990 

2.83e-02 

0.995 

optimal 

uy 

2.32e-01 

1.17e-01 

0.988 

5.85e-02 

1.000 

2.93e-02 

0.998 

optimal 

3 

u 

1.72e-04 

l.lle-05 

3.954 

6.99e-07 

3.989 

4.37e-08 

4.000 

optimal 

u X 

4.12e-03 

5.00e-04 

3.043 

6.17e-05 

3.019 

7.69e-06 

3.004 

optimal 

uy 

4.46e-03 

5.42e-04 

3.041 

6.71e-05 

3.014 

8.35e-06 

3.006 

optimal 

Tables  10  and  11  show  the  convergence  rates  of  pi  and  p3  approximations  on  rectangular  and  triangular 
meshes,  respectively.  Not  surprisingly,  the  convergence  rates  are  optimal  for  all  cases. 

In  Fig.  7,  the  solution  u  and  its  gradient  ux  along  the  horizontal  line  y  =  0  are  shown  together  with  the 
analytical  solutions  at  t  =  0.05.  To  visually  compare  the  accuracy  between  pi  and  p3  results,  we  intentionally 
choose  a  coarse  20  x  20  mesh.  As  can  be  seen,  the  pi  solution  is  not  accurate  enough  to  resolve  the  local 
extrema.  By  contrast,  the  p3  solution  lies  on  top  of  the  analytical  solution. 


4.2.2  2-D  Viscous  Burgers  Equation 


Finally,  the  following  2-D  viscous  Burgers  equation  is  solved  using  DG-CVS. 


du  du  du  ( d2u 

di+t‘di+Udy-‘'W 


0  <  x,  y  <  25 


with  the  analytical  solution  given  as 


2 

^exact  —  "  x-xc+y-yc-2t 

e  -  + 1 


(32) 


where  (xc,yc)  =  (0,0)  is  a  constant  location. 

The  1-D  version  of  this  case  was  presented  in  [18].  This  case  is  constructed  such  that  the  original  wave 
is  propagated  without  changing  shape  under  the  effect  of  both  nonlinear  advection  and  linear  diffusion.  The 
initial  solution  at  t  =  0  and  the  boundary  conditions  on  all  four  boundaries  are  provided  by  the  analytical 
solution.  Therefore  the  boundary  conditions  are  time  dependent. 

We  first  conduct  the  grid  convergence  study  on  this  nonlinear  advection- diffusion  case.  In  the  study, 
v  =  2.5  is  chosen  and  the  simulation  is  run  up  to  t  =  10.  In  the  current  study,  the  time  step  is  chosen  as 
St  =  crmin(/i/a,  h2/v)  where  h  is  the  local  mesh  size  and  a  =  0.2  for  the  pi  case  and  a  =  0.1  for  the  p3 
case.  For  the  purpose  of  determining  the  convergence  rates,  this  choice  of  time  steps  may  not  be  appropriate 
since  advection  and  diffusion  have  different  time  scales.  A  more  appropriate  approach  may  be  the  operator 
splitting  method  where  the  advection  and  the  diffusion  are  treated  with  different  time  steps.  Tables  12 
and  13  show  the  convergence  rates  of  pi  and  p3  on  rectangular  and  triangular  meshes,  respectively.  The 
convergence  rates,  though  not  as  neat  as  those  for  the  pure  diffusion  equation,  are  still  close  to  optimal. 

To  further  demonstrate  the  accuracy  of  DG-CVS  for  various  values  of  z/,  the  following  three  cases  are 
simulated  on  both  qua-50  and  tri-50  meshes  where  the  size  of  the  mesh  is  indicated  by  5x  =  0.5: 
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Solution  along  line  y=0,  v=1.0,  T=0.05 


Solution  along  line  y=0,  v=1.0,  T=0.05 


Solution  along  line  y=0,  v=1.0,  T=0.05 


Solution  along  line  y=0,  v=1.0,  T=0.05 


Figure  7:  Comparison  between  pi  and  p3  solutions  on  the  20  x  20  mesh  at  y  =  0  and  t  =  0.05  for  the  case 
of  the  fundamental  solution  of  the  2-D  heat  equation. 


Table  12:  Solution  of  2-D  viscous  Burgers  equation  on  rectangular  meshes,  v  =  2.5,  t  =  10.  Numerical 
convergence  order  determined  by  L 2  norm. 


p 

variable 

qua- 10 

qua- 20 

qua-40 

qua- 80 

L2(e) 

L2(e) 

order 

L2(e) 

order 

i2(e) 

order 

comments 

1 

u 

ux  and  uy 

1.38e-02 

1.74e-02 

3.15e-03 

8.52e-03 

2.131 

1.030 

7.71e-04 

4.21e-03 

2.031 

1.017 

1.92e-04 

2.09e-03 

2.006 

1.010 

optimal 

optimal 

3 

u 

ux  and  uy 

2.83e-04 

7.11e-04 

1.69e-05 

9.12e-05 

4.066 

2.963 

l.lle-06 

1.19e-05 

3.928 

2.938 

7.88e-08 

1.61e-06 

3.816 

2.886 

optimal 

optimal 

Table  13:  Solution  of  2-D  viscous  Burgers  equation  on  triangular  meshes,  v  =  2.5,  t  =  10.  Numerical 
convergence  order  determined  by  L2  norm. 


tri-10 

tri-20 

tri-40 

tri-80 

P 

variable 

L2(e) 

L2(e) 

order 

L2(e) 

order 

L2(e) 

order 

comments 

1 

u 

8.08e-03 

2.19e-03 

1.883 

5.49e-04 

1.996 

1.38e-04 

1.992 

optimal 

ux 

1.54e-02 

7.79e-03 

0.983 

3.86e-03 

1.013 

1.93e-03 

1.000 

optimal 

uy 

1.54e-02 

7.79e-03 

0.983 

3.86e-03 

1.013 

1.93e-03 

1.000 

optimal 

3 

u 

1.70e-04 

1.21e-05 

3.812 

1.00e-06 

3.597 

1.03e-07 

3.279 

optimal 

ux 

4.16e-04 

5.31e-05 

2.970 

6.44e-06 

3.044 

7.96e-07 

3.016 

optimal 

uy 

4.16e-04 

5.34e-05 

2.962 

6.48e-06 

3.043 

8.02e-07 

3.014 

optimal 
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Solution  along  line  x-y=0,  v=0.05,  T=20 


Solution  along  line  x-y=0,  v=0.5,  T=20 


Solution  along  line  x-y=0,  v=2.5,  T-20 
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0  4  8  12  16  20  24  28  32  36 


V2x 


V2x 


V2x 


Figure  8:  pi  solution  of  2-D  viscous  Burgers  equation  along  the  diagonal  line  x  —  p  =  0  at  £  =  20.  Left: 
v  =  0.05,  middle:  v  =  1.0,  and  right:  v  —  2.5. 


Solution  along  line  x-y=0,  v=0.5,  T=20 


Solution  along  line  x-y=0,  v=0.5,  T=20 


V2x 


V2x 


Figure  9:  ux  or  uy  of  the  2-D  viscous  Burgers  equation  with  v  =  0.5  along  the  diagonal  line  at  t  =  20.  Left: 
pi  solution  and  right:  p3  solution. 


•  short  wave  ( 5x/v  =  10,  i.e.  v  =  0.05). 

•  medium  wave  (Sx/v  =  1,  i.e.  v  =  0.5). 

•  long  wave  ( dx/v  =  0.2,  i.e.  v  =  2.5). 

Figure  8  plots  the  pi  solution  along  the  diagonal  line  x  —  y  =  0  together  with  the  exact  solution.  As  can 
be  seen,  all  types  of  waves  have  been  captured  accurately  in  terms  of  both  location  and  magnitude.  In  the 
case  of  the  short  wave,  no  limiter  is  applied  and  therefore  slight  overshoot  and  undershoot  can  be  seen. 
The  p3  solution  is  not  shown  since  no  much  visual  difference  can  be  seen  between  the  p3  solution  and  the 
pi  solution.  Figure  9  compares  the  computed  solution  gradients  ( ux  or  uy)  with  the  exact  solution  for  the 
case  of  v  =  0.5.  As  can  be  seen,  the  p3  solution  is  superior  to  the  pi  solution  in  resolving  the  local  sharp 
extremum  of  the  solution  gradient. 

5  Major  Progress  #2:  Solving  Euler  Equations  Involving  Curved 
Boundaries 

Progress  has  been  made  to  solve  flow  problems  involving  curved  boundaries.  Here  two  examples  are  presented 
to  demonstrate  the  capability  of  handling  curved  boundaries. 

5.1  Flow  around  a  Circular  Cylinder 

The  first  example  is  a  subsonic  flow  with  M  =  0.1  around  a  circular  cylinder.  The  far-held  boundary  is 
located  at  30  diameters  from  the  center  of  the  circle.  The  mesh  contains  24  x  10  quadrilateral  cells.  Figure 
10  shows  the  Mach  number  distributions  from  the  second-order  (pi),  third  order  (p2)  and  fourth  order  (p3) 
simulations,  respectively.  Obviously,  high-order  solutions  are  superior  to  lower  order  ones  on  the  same  coarse 
mesh.  Figure  11  compares  the  steady-state  pressure  distribution  on  the  cylinder  surface  with  the  analytical 
solution.  Again,  p3  solution  is  the  most  accurate. 
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Figure  10:  Mach  number  field  of  the  subsonic  flow  (M  =  0.1)  around  a  circular  cylinder  on  a  24  x  10 
quadrilateral  mesh.  Left:  pi  solution,  middle:  p2  solution,  and  right:  p3  solution. 


pressure  on  the  cylinder  surface  (M=0.1) 


pressure  on  the  cylinder  surface  (M=0.1) 


pressure  on  the  cylinder  surface  (M=0.1) 


x 


x 


x 


Figure  11:  Pressure  on  the  surface  of  the  circular  cylinder  on  a  24  x  10  quadrilateral  mesh.  Left:  pi  solution, 
middle:  p2  solution,  and  right:  p3  solution. 
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Figure  12:  Subsonic  flow  around  NACA0012  airfoil.  Left:  pressure  field.  Right:  Mach  number  field. 

5.2  Flow  around  NACA0012  Airfoil 

The  second  example  is  a  subsonic  flow  with  Mach  0.63,  AoA  =  2°  around  the  NACA0012  airfoil  with  curved 
boundary.  The  mesh  contains  1536  quadrilateral  cells  with  64  cells  on  the  airfoil  surface.  Figure  12  shows 
the  p3  solution.  As  can  be  seen,  the  contour  lines  are  smooth  near  the  boundary  including  the  trailing  edge. 


6  Major  Progress  #3:  Solving  Shallow  Water  Equations 


The  shallow  water  equations  (SWE)  are  widely  used  as  the  mathematical  model  to  numerically  simulate 
the  dam  break,  river  inundation,  failure  of  levees  and  tide  of  ocean  in  coastal  and  civil  engineering.  When 
the  characteristic  vertical  velocity  is  small  in  comparison  with  the  characteristic  horizontal  velocity,  which 
happens  when  the  characteristic  vertical  length  scale  is  much  smaller  than  the  characteristic  horizontal  length 
scale,  the  incompressible  Navier-Stokes  equations  can  be  simplified  to  the  shallow  water  equations  (SWE) 
by  depth  averaging  the  NSE  in  the  vertical  direction. 

The  frictionless  shallow  water  equation  can  be  expressed  in  the  following  conservative  form 


du  df  dg  , 

m+d^  +  d^  =  r{u) 

in  which  the  conservative  vector,  flux  vectors  and  source  vector  are 


H 

Hu 

Hv 

0 

u  = 

Hu 

Hv 

,f  = 

Hu 2  +  gH2/ 2 
Huv 

,  g  = 

Huv 

Hv2  +  gH2/ 2 

,  and  r  = 

qH  Sqx 

gHSoy 

(33) 


respectively.  Here,  the  total  water  depth  iL(x,t)  =  £(x,  £)  +  d(x)  where  £(x,  £)  is  the  free-surface  elevation 
and  d(x)  is  the  still  water  depth,  u  and  v  are  the  x-  and  //-component  of  the  depth- averaged  velocity,  g  is 
the  acceleration  due  to  gravity.  Sqx  and  Soy  are  the  bottom  slopes  in  x-  and  //-directions,  respectively. 

Equation  (33)  is  a  hyperbolic  nonlinear  system.  There  is  close  mathematical  and  physical  analogy  be¬ 
tween  the  shallow  water  flows  and  compressible  flows.  The  hydraulic  jumps  and  bores  are  analogous  to  the 
stationary  and  moving  shock  waves  in  compressible  gas  flows.  Therefore,  the  numerical  methods  used  to 
solve  the  SWE  often  mimic  those  for  solving  the  compressible  Euler  equations. 

The  current  DG-CVS  method  provides  a  Riemann-solver-free  alternative  approach  to  solve  the  shallow 
water  equations.  Here,  several  numerical  examples  are  provided  to  test  the  accuracy  of  the  DG-CVS  based 
solver  on  solving  the  shallow  water  equations.  In  all  the  simulations,  the  reference  velocity  is  chosen  according 
to 

^ref  =  ref 
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Figure  13:  Water  elevation  carpet  view  of  the  2-d  shocktube  problem.  Left:  t  =  0.0,  middle:  t  =  1.0  and 
right:  t  =  2.0. 


where  g  is  the  acceleration  due  to  gravity  and  Lref  is  some  reference  length.  Nondimensionalization  based 
on  these  quantities  leads  to  that  the  nondimensional  g  is  unity.  In  addition,  all  simulations  assume  flat  water 
bed. 

6.1  2-D  Dam  Break  Problem 

The  first  case  is  a  two-dimensional  dam  break  problem  which  is  analogous  to  the  shock  tube  problem  used 
to  verify  a  compressible  Euler  solver.  The  computational  domain  is  a  [—5,5]  x  [—5,5]  rectangular  region. 
At  t  =  0,  the  water  height  is  3  within  the  lower  left  corner  [—5,0]  x  [—5,0],  and  the  rest  of  the  domain  is 
filled  with  water  of  height  1.  The  left  and  the  bottom  boundaries  are  symmetry  ones.  Figure  13  shows  the 
carpet  view  of  the  water  elevation  at  three  time  instants.  The  solution  is  obtained  using  the  fourth  order 
basis  functions  on  the  99  x  99  mesh.  As  can  be  seen,  the  complex  wave  structures  have  been  captured  in 
the  simulation.  Note  that  no  any  type  of  solution  limiting  is  employed  in  the  simulation.  Therefore,  slight 
wiggles  can  be  seen  behind  the  water  jump.  Figure  14  shows  the  water  elevation  and  momentum  along  one 
of  the  symmetry  boundaries  at  times  up  to  t  =  2.0.  The  water  jump  discontinuity  is  captured  sharply. 

6.2  Gaussian  Pulse  Problem 

The  second  case  is  to  simulate  the  evolution  of  a  Gaussian  water  depth  perturbance.  This  problem  is  indeed 
one-dimensional  but  is  simulated  on  a  2-D  computational  domain  as  shown  in  Fig.  15.  The  fourth  order  (p3) 
simulation  is  performed  again.  Figures  15  and  16  shows  similar  information  to  those  in  Figs.  13  and  14  in 
the  previous  example.  As  can  be  seen,  the  two  waves  propagate  to  two  opposite  directions  and  get  steeper 
and  steeper. 

6.3  2-D  Dam  Break  Problem 

The  last  example  is  a  more  realistic  classical  benchmark  problem  to  verify  a  shallow  water  equation  solver. 
Figure  17  depicts  the  domain  geometry.  Initially,  two  reservoirs  are  separated  by  a  lock  gate  which  is  of  75 
meters  wide.  The  water  levels  are  10  meters  and  5  meters,  respectively.  The  reservoir  on  the  left  has  higher 
water  level.  At  t  =  0,  the  lock  gate  is  opened.  The  evolution  of  the  water  level  is  simulated  using  the  second 
order  (pi)  DG-CVS  solver.  The  resolution  of  the  computational  mesh  is  2.5  meters.  Figure  18  shows  the 
contour  lines  and  carpet  view  of  the  water  level  at  t  =  7.2  seconds. 
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H  at  t  =  0 


Figure  14:  Water  element  and  momentum  along  the 
and  right:  momentum. 


Hu  at  t  =  0 


Figure  15:  Carpet  view  of  the  evolution  of  Gaussian  pulse  depth  perturbation.  Clockwisely,  t  =  0,  t  =  0.5, 
t  =  1.0  and  t  =  2.0. 


7  Major  Progress  ^4:  Solving  Level  Set  Equation 

The  DG-CVS  solver  has  been  also  extended  to  solve  the  level-set  equation  to  accurately  resolve  the  interface 
between  fluids.  In  the  level  set  method,  the  interface  is  represented  by  the  zero  level  set  C{t)  =  {x| ^(x,  t)  =  0} 
of  a  level  set  function  ^(xT).  'ip  is  often  taken  to  be  the  signed  distance  function  to  the  interface.  The  level 
set  equation  can  be  expressed  as 


'ipt  +  v  •  =  0  (34) 

where  v  is  the  velocity  field  which  evolves  the  interface.  Equation  (34)  can  be  rewritten  in  conservative  form 

ipt  +  V  •  (V>v)  =  'ipV  •  v  (35) 


Several  examples  are  provided  here  to  demonstrate  the  capability  of  the  DG-CVS  solver  to  solve  the 
level-set  equation.  In  all  the  examples,  the  initial  signed  distance  field  is  narrowed  by  the  following 


ip 


— e  if  ip  <  —e 
<  ip  if  —  e  <  'ip  <  e 
e  if  'ip  >  e 


wheree  is  chosen  to  be  0.1.  This  strategy  is  based  on  the  assumption  that  the  interface  evolution  often  occurs 
within  a  narrow  band. 


7.1  Rotation  of  an  off-center  circular  fluid  body 

The  first  example  is  about  the  rotation  of  a  circular  fluid  body.  The  circular  fluid  body  of  radius  0.15  is 
centered  at  the  location  (0.5,  0.75).  The  circle  is  rotated  by  a  velocity  field  v  =  (u,v)  given  by: 

u  =  2tt(i/  —  0.5),  and  v  =  —  2ir(x  —  0.5).  (36) 

Such  a  velocity  field  completes  one  revolution  per  unit  time.  Figure  19  is  the  initial  signed  distance  field 

where  the  zero  level  set  is  also  indicated.  The  simulation  is  done  on  a  40  x  40  mesh.  Figures  20  and  21 

compares  the  pi  and  p3  solutions  at  various  times.  Again,  p3  solution  preserves  the  circular  shape  better 
than  the  pi  solution. 
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Figure  16:  Line  plots  of  the  evolution  of  Gaussian  pulse  depth  perturbation. 
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Figure  18:  Contour  lines  and  the  carpet  view  of  the  water  elevation  of  the  2-D  dam  break  problem  at  t  =  7.2s. 


Figure  19:  Initial  signed  distance  field  of  an  off-center  circular  fluid  body. 
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(c)  t  =  2.0. 


(d)  t  =  3.0. 


Figure  20:  Rotation  of  an  off-center  circular  fluid  brjl'ly  on  a  40  x  40  rectangular  mesh.  Left  column:  pi 
solutions  and  right  column:  p3  solutions,  (to  be  continued) 


(a)  t  =  4.0. 


(b)  t  =  5.0. 


Figure  21:  Rotation  of  an  off-center  circular  fluid  body  on  a  40  x  40  rectangular  mesh.  Left  column:  pi 
solutions  and  right  column:  p3  solutions,  (continued) 
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Figure  22:  Zalesak  disk.  Initial  level  set  is  the  signed  distance  to  the  disk. 


7.2  Rotation  of  the  slotted  Zalesak  disk 

This  another  well  known  example  is  the  rigid  rotation  of  the  Zalesak  disk  [19].  The  initial  interface  is  a  circle 
centered  at  (0.5,  0.75)  with  a  notch  of  width  and  depth  0.05  and  0.25,  respectively.  The  velocity  field  is  the 
same  as  the  previous  case,  i.e.  (36). 

The  initial  level  set  function  is  initialized  as  the  signed  distance  to  the  disk  (Fig.  22).  The  zero  level 
set  indicates  the  interface  location.  Figure  23  shows  the  fourth  order  (p3)  solution  using  the  solver  based 
on  DG-CVS.  The  computational  mesh  contains  50  x  50  rectangular  cells.  The  time  step  in  the  simulation 
is  St  =  0.00125.  800  time  steps  were  finish  one  revolution.  Compared  to  the  exact  solution,  the  smearing  is 
very  small  which  indicates  the  small  numerical  dissipation  of  this  high  order  scheme. 


7.3  Shearing  of  a  circular  bubble 

The  last  example  is  to  simulate  the  deformation  of  a  circular  bubble  of  radius  0.15  centered  at  (0.5,  0.75). 
The  circle  is  deformed  due  to  the  following  velocity  field  [2C] : 

u  =  —  sin(27 ry)  sin2  (to)  cos(7 rt/T),  and  v  =  sin(2TO)  sin2  (7 xy)  cos(7r t/T)  (37) 


which  are  obtained  from  (— <9T/ch/,  d^/dx)  where  T  is  the  following  stream  function 


T  =  —  sin2 (to)  sin2 (7 xy). 

IX 


At  t  =  T/2,  the  velocity  starts  to  reverse  direction  and  at  time  t  =  T,  the  original  circle  is  expected  to  be 
restored. 

The  velocity  field  given  by  (37)  represents  a  vortex  which  stretch  the  circular  fluid  body  into  a  thin 
filament  toward  the  vortex  center.  The  simulation  is  performed  on  a  relatively  coarse  40  x  40  mesh.  Here 
T  =  4.  Figures  24  and  25  shows  the  fourth  order  (p3)  evolution  of  the  signed  distance  field  and  the  zero 
level  set  at  various  times.  At  t  —  0.5 T  =  4,  the  bubble  starts  to  restore  its  shape.  At  t  m  4.0,  the  bubble 
recovers  its  original  circular  shape.  The  comparison  with  the  exact  solution  shows  the  accuracy  of  the  current 
method. 


8  Major  Progress  7^5:  Solving  Moving  Mesh  Problems 

The  ultimate  goal  of  developing  DG-CVS  based  solvers  is  to  solve  fluid  flow  problems  involving  moving 
boundaries.  Due  to  the  space-time  formulation  of  DG-CVS,  DG-CVS  automatically  satisfies  the  so-called 
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(a)  t  =  0.25.  (b)  t  =  0.5. 


(c)  t  =  1.0.  (d)  t  =  1.5. 


(e)  t  =  2.0.  (f)  t  =  2.5. 

Figure  24:  p3  solution  of  the  shearing  of  a  circular  bubble  by  a  vortex  on  a  40  x  40  mesh.  T  =  4. 
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(a)  t  =  3.0.  (b)  t  =  3.5. 


Figure  25:  p3  solution  of  the  shearing  of  a  circular  bubble  by  a  vortex  on  a  40  x  40  mesh.  T  =  4. (continued) 
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Figure  26:  Meshes  used  in  the  simulation.  Left:  20  x  20  rectangular  mesh  qua- 20.  Right:  unstructured 
triangular  mesh  tri-20  with  similar  edge  resolution. 


geometric  conservation  law  (GCL)  by  incorporating  the  mesh  motion  effect  into  the  unit  normal  of  the 
space-time  side  surfaces  of  the  conservation  element  during  the  enforcement  of  space-time  fluxes.  Therefore, 
any  extra  sophisticated  techniques  used  to  enforce  the  GCL  in  other  semi-discrete  high-order  methods  can 
be  avoided. 

A  simple  test  is  presented  here  to  demonstrate  the  moving- mesh  capability  of  the  current  DG-CVS 
solver.  This  test  is  an  isentropic  vortex  advection  problem  on  a  2D  square  domain  [0, 10]  x  [0, 10].  The  initial 
conditions  are  given  as  an  isentropic  vortex  with  the  center  at  (5,  5),  i.e., 

u{x,  y,  0)  =  1  —  J_e0-5(1_r2)  (y  -  5)  (38a) 

v(x,y,  0)  =  1  +  ;^-e°'5(1-r  —  5)  (38b) 

27 T 

T(x,  y,  0)  =  1  —  A-T^e1_r2  (38c) 

S(x,y,0)*l  (38d) 

where  u,  u,  T,  and  S  are  x- velocity,  2/- velocity,  temperature  and  entropy,  respectively,  e  =  5  represents  the 
vortex  strength  and  r2  =  (x  —  5) 2  +  (y  —  5)2.  The  periodic  boundary  conditions  on  both  directions  are 
assumed.  The  density  p  and  the  pressure  P  can  be  obtained  via 

P(x,y,  0)  =  (  >  P(x,y,0)  =  p{x,y,0)T(x,y,0) 

\S{x,y,0)J 

where  7  =  1.4  is  the  ratio  of  specific  heats.  It  can  be  shown  that  the  Euler  equations  with  the  above  initial 
and  boundary  conditions  allows  an  exact  solution  which  is  the  initial  solution  advected  with  the  speed  (1,1) 
in  the  diagonal  direction.  With  periodic  boundary  conditions,  the  vortex  will  come  back  to  the  center  of  the 
domain,  its  original  location,  at  t  =  10. 

Both  rectangular  and  triangular  meshes  are  tested.  Since  we  are  going  to  present  the  fourth  order  (p3) 
solutions  in  this  paper,  relatively  coarse  meshes  are  used.  The  meshes  are  designated  as  rec-20  and  tri-20 
(cf.  Fig.  26)  where  each  boundary  of  the  domain  contains  equally  spaced  20  edges.  The  interior  elements  of 
the  triangular  meshes  are  unstructured  and  generated  by  EasyMesh  [  ]. 

In  the  simulation,  the  grid  points  on  the  boundaries  are  fixed  and  the  interior  grid  points  move  according 
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Figure  27:  Comparison  of  the  p3  density  field  at  t  =  8  of  the  vortex  advection  problem  on  mesh  qua-20. 
Left:  stationary  mesh.  Right:  moving  mesh. 


to 


x(£,  77,  t)  =  £  +  sin(0.27r£)  sin(0.27r?7)  sin(— 0.27r£)  (39a) 

t/(£,  77,  t)  =  77  +  sin(0.27r^)  sin(0.27r?7)  sin(— 0.27r£)  (39b) 

where  (£,77)  is  the  spatial  coordinate  in  the  fixed  reference  domain  which  is  the  same  as  the  physical  domain 
at  t  =  0.  The  above  mesh  motion  equations  also  ensure  that  the  grid  points  on  the  vertical  and  horizontal 
center  lines  (£  =  5  and  77  =  5,  respectively)  are  not  moving  with  time. 

For  fair  comparison,  the  time  step  used  in  all  simulations  is  taken  as  St  =  0.02,  regardless  of  whether  the 
mesh  is  rectangular  or  triangular  and  whether  the  mesh  is  stationary  or  moving.  Only  the  fourth  order  (p3) 
results  will  be  provided. 

Figure  27  shows  the  comparison  of  the  density  fields  at  t  =  8  on  the  stationary  and  the  moving  quadrilat¬ 
eral  meshes.  As  can  be  seen,  the  solution  on  the  moving  mesh  is  visually  as  smooth  as  the  one  on  stationary 
mesh.  Figure  28  shows  similar  density  fields  on  triangular  meshes. 

Figures  29  and  30  compare  the  density  distribution  at  t  =  10  along  the  horizontal  center  line  (y  =  5)  of 
the  domain  with  the  exact  solution  for  all  simulations.  The  accuracy  is  clearly  demonstrated. 

Finally,  we  present  the  quantitative  error  comparison  between  all  simulations  in  Table  14.  Here,  A,  I2 
and  Zqo  errors  are  evaluated  at  the  discrete  vertices  of  the  mesh.  The  L2  error  is  the  integrated  error  over 
the  top  surfaces  of  the  vertex-level  CE.  Recall  that  all  simulations  use  the  same  time  step  St  =  0.02.  As 
can  be  seen,  for  the  same  qua-20  or  tri-20  meshes,  the  errors  on  moving  meshes  are  larger  than  those  on 
stationary  meshes.  The  is  expected  because  the  current  prescribed  mesh  motion  is  independent  from  the 
solution.  Mesh  motion  causes  the  loss  of  grid  resolution  compared  with  the  uniform  stationary  mesh  (cf. 
Figs.  27  and  28).  In  addition,  the  errors  on  triangular  meshes  are  smaller  than  those  on  rectangular  meshes. 
This  is  because  there  are  more  cells  on  the  triangular  mesh  even  if  its  edge  resolution  is  comparable  to  that 
of  the  rectangular  mesh.  Therefore,  triangular  meshes  yield  more  accurate  results  at  the  price  of  higher 
computational  cost. 
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Figure  28:  Comparison  of  the  p3  density  field  at  t  =  8  of  the  vortex  advection  problem  on  mesh  tri-20. 
Left:  stationary  mesh.  Right:  moving  mesh. 


t=10,  rectangular  stationary  mesh 


t=10,  rectangular  moving  mesh 


Figure  29:  Comparison  of  the  p3  density  distribution  of  the  vortex  advection  problem  at  t  =  10  along  the 
horizontal  center  line  of  the  domain  on  mesh  qua- 20.  Left:  stationary  mesh.  Right:  moving  mesh. 


t=10,  triangular  stationary  mesh  t=10,  triangular  moving  mesh 


Figure  30:  Comparison  of  the  p3  density  distribution  of  the  vortex  advection  problem  at  t  =  10  along  the 
horizontal  center  line  of  the  domain  on  mesh  tri-20.  Left:  stationary  mesh.  Right:  moving  mesh. 
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Table  14:  Comparison  of  the  p3  density  errors  of  the  isentropic  vortex  advection  problem  on  fixed  and  moving 
meshes. 


fixed  qua-20  mesh 

moving  qua-20  mesh 

fixed  tri-20  mesh 

moving  tri-20  mesh 

l\ -error 

4 . 32e-05 

2 . 54e-04 

2 . 13e-05 

8 . 41e-05 

I2  -error 

1 . 19e-04 

4 . 22e-04 

3 . 77e-05 

1 . 80e-04 

Zoo -error 

1 . 50e-03 

3 . 36e-03 

4 . 82e-04 

1 . 68e-03 

L2 -error 

1 . 12e-04 

4 . 58e-04 

4 . 15e-05 

1 . 84e-04 
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